% =========================================================================
% Schniepp Lab, 2018-2021
% The second part to carry out the unpolarized Raman spectra fitting
% Steps:
% 3. Model polypeptide Raman spectra interpolation and normalization
% 4. Calculate the unpolarized Raman spectrum
% 5. Fitting
% 6. Result plotting
% =========================================================================

% Initialize the MATLAB workspace
clear;
load('./Data/Cleaned Polypeptide Data.mat');

% =========================================================================
% 3. Model polypeptide Raman spectra interpolation and normalization
% =========================================================================

% Create the x values used for interpolation
xR1 = [920:1:1800]';

% Raman spectra of model poly peptides
% alpha - Poly-alanine
Unpolar_Raman_aPA_Fanconi1973_G(:,1) = xR1;
Unpolar_Raman_aPA_Fanconi1973_G(:,2) = interp1(Unpolar_Raman_aPA_Fanconi1973_P(:,1),Unpolar_Raman_aPA_Fanconi1973_P(:,2),xR1);
Unpolar_Raman_aPA_Fanconi1973_G(isnan(Unpolar_Raman_aPA_Fanconi1973_G)) = 0;

% beta - Poly-alanine
Unpolar_Raman_bPA_G(:,1) = xR1;
Unpolar_Raman_bPA_G(:,2) = interp1(Unpolar_Raman_bPA_P(:,1),Unpolar_Raman_bPA_P(:,2),xR1);
Unpolar_Raman_bPA_G(isnan(Unpolar_Raman_bPA_G)) = 0;

% Glycine - I
Unpolar_AmideI_Raman_GlyI_G(:,1) = xR1;
Unpolar_AmideI_Raman_GlyI_G(:,2) = interp1(Unpolar_AmideI_Raman_GlyI_P(:,1),Unpolar_AmideI_Raman_GlyI_P(:,2),xR1);
Unpolar_AmideI_Raman_GlyI_G(isnan(Unpolar_AmideI_Raman_GlyI_G)) = 0;

% Glycine - II
Unpolar_AmideI_Raman_GlyII_G(:,1) = xR1;
Unpolar_AmideI_Raman_GlyII_G(:,2) = interp1(Unpolar_AmideI_Raman_GlyII_P(:,1),Unpolar_AmideI_Raman_GlyII_P(:,2),xR1);
Unpolar_AmideI_Raman_GlyII_G(isnan(Unpolar_AmideI_Raman_GlyII_G)) = 0;

% Rough magnitude alignment
aPA = Unpolar_Raman_aPA_Fanconi1973_G(:,2)*10;
bPA = Unpolar_Raman_bPA_G(:,2)*10;
gIP = Unpolar_AmideI_Raman_GlyI_G(:,2)*10;
gIIP = Unpolar_AmideI_Raman_GlyII_G(:,2)*10;

% Crop the polypeptide data to leave the range of interest
xR1 = xR1(232:785);
aPA = aPA(232:785);
bPA = bPA(232:785);
gIP = gIP(232:785);
gIIP = gIIP(232:785);

% Consolidate them into a single variable
pep_data(:,1) = aPA;
pep_data(:,2) = bPA;
pep_data(:,3) = gIP;
pep_data(:,4) = gIIP;

% Calculate the normalization factor of the polypeptide spectra to make 
% them have the same area underneath the curve
target = 10^5;
a1 = trapz(pep_data(:, 1));
a2 = trapz(pep_data(:, 2));
a3 = trapz(pep_data(:, 3));
a4 = trapz(pep_data(:, 4));

r1 = target / a1;
r2 = target / a2;
r3 = target / a3;
r4 = target / a4;

% Normalize the polypeptide using the normalization factor calculated above
pep_data(:, 1) = pep_data(:, 1) * r1;
pep_data(:, 2) = pep_data(:, 2) * r2;
pep_data(:, 3) = pep_data(:, 3) * r3;
pep_data(:, 4) = pep_data(:, 4) * r4;

% =========================================================================
% 4. Calculate the unpolarized Raman spectrum
% =========================================================================

% Load the experimental Raman polarized data
load('./Data/Raman_Exp_Data.mat');

% Interpolation for experimental data
new_freq = [1151:1:1704]';
new_xx = interp1(freq, xx, new_freq);
new_yy = interp1(freq, yy, new_freq);
new_zz = interp1(freq, zz, new_freq);
new_zx = interp1(freq, zx, new_freq);
new_xz = interp1(freq, xz, new_freq);

% Assume xy = yx = (xx + yy) / 2
% Assume yz = zy = xz / 1.5 = zx / 1.5 = (xz + zx) / 3
% Unpolarized spectrum = xx + yy + zz + xz + zx + yz + zy + xy + yx
% = xx + yy + zz + xz + zx + 2 * (xz + zx) / 3 + (xx + yy)
% = 2 * xx + 2 * yy + zz + 5 * xz / 3 + 5 * zx / 3;
exp_data = 2 * new_xx + 2 * new_yy + new_zz + 5 * new_zx / 3 + 5 * new_xz / 3;
spectra = 'all';

% Exclude the functional group peaks
indices = [1:360, 481:554];
Range_string = 'NoFuncGroup';
pep_data_cut = pep_data(indices, :);
exp_data_cut = exp_data(indices, :);

% =========================================================================
% 5. Fitting
% =========================================================================

% Initialize the coefficents for each polypeptide component and define the
% fitting lower and upper bounds
% [aPA, bPA, gIP, gIIP]
coeff_guess = [1.0, 1.0, 1.0, 1.0];
lb = [0, 0, 0, 0];
ub = [10, 10, 10, 10];

pep_comp = 'bPA+gIP+gIIP+aPA';

% Carry out the actual fitting
[coeff, resnorm, resid, J] = optimizeCoeffwithMore(coeff_guess, pep_data_cut, exp_data_cut, lb, ub, 'Cal', 1000000);
ci = nlparci(coeff, resid, 'jacobian', J);
error = norm(resid);

% =========================================================================
% 6. Result plotting
% =========================================================================

% Plot the calculated unpolarized Raman spectrum
figure();
plot(new_freq, exp_data, 'b', 'linewidth', 1.5);
hold on;
grid on;
set(gca,'xlim',[1100 1800]);
set(gcf,'position',[180, 39, 1145, 746]);
set(gca,'fontsize',15);

% Plot the polypeptide fitting result
n_sR = pep_data(:, 1) * coeff(1) + pep_data(:, 2) * coeff(2) + pep_data(:, 3) * coeff(3) + pep_data(:, 4) * coeff(4);
plot(new_freq, n_sR, 'r');
plot(new_freq, pep_data(:, 1) * coeff(1), 'k');
plot(new_freq, pep_data(:, 2) * coeff(2), 'm');
plot(new_freq, pep_data(:, 3) * coeff(3), 'g');
plot(new_freq, pep_data(:, 4) * coeff(4), 'c');
t_string = strcat('Raman,', spectra, ',', pep_comp, ',', Range_string);
title(t_string);
set(gca, 'xlim', [1150 1704]);
s_aPA = strcat('aPA,', num2str(coeff(1),3));
s_bPA = strcat('bPA,', num2str(coeff(2),3));
s_gIP = strcat('gIP,', num2str(coeff(3),3));
s_gIIP = strcat('gIIP,', num2str(coeff(4),3));
legend('Exp','Fit',s_aPA,s_bPA,s_gIP,s_gIIP, 'Location', 'northwest');